%pylab inline
import matplotlib.pyplot as plt
from ipywidgets import *  # az interaktivitásért felelős csomag

import collections

Populating the interactive namespace from numpy and matplotlib

Calculation of the Landau levels using tight-binding approximation

for a finite width of strip in magnetic field perpendicular to the plane of the square lattice.

Landau gauge: the vektor potential $\mathbf{A}$ is along the $x$ axis, but proportional to the coordinate $y$, ie, $\mathbf{A} = (-B y, 0, 0)$.

Using the Peierls substitution and in Landau gauge the system is translation invariant along the $x$ axis, then $\left[H,p_x\right] = 0$

Parameter dim -- > is the width of a slice, the number of atoms in the $y$ direction.

The code calculates the Landau levels.

Hamilton operator:

$H = H_0 + H_1 e^{i k} + H_1^+ e^{-i k}$,

where $H_0$ is the Hamiltonian for one slice and $H_1$ is the connection of one slice with the next slice to the right.

Current operator (depending on $y$ but along the $x$ direction):

$\hat{j}_x = \frac{\partial H}{\partial k} = i k \left( H_1 e^{i k}- H_1^+ e^{-i k}\right)$

$j_x(y) = \psi^+_n(y) \frac{\partial H}{\partial k} \psi_n(y) $

# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('Fig_square_lattice_srtip_Landau_edge_state.pdf');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 8  #    12 nagy abrahoz
yfig_meret= 6    #   12 nagy abrahoz
xyticks_meret= 20  #  20 nagy abrahoz
xylabel_meret= 20  #  30 nagy abrahoz
legend_meret= 20   #  30 nagy abrahoz

def Htot(k, theta, dim):
    
    # dim ---> y iranyban az elemi cellak szama
    
    H0 = matrix(diag((-1)*ones(dim-1),1),complex)
    H0 = H0 + H0.H
    
    H1 = matrix(zeros((dim,dim),complex))
    for i in range(dim):
        #H1[i,i] = -exp(1j*i*theta)
        H1[i,i] = -exp(-1j*(i-dim/2)*theta)  
        # csak az van, hogy a diszperzios relacio k-ban eltolodik
        
    Htot = H0 + exp(1j*k)*H1 + exp(-1j*k)*H1.H
    Hcurrent = 1j*k*(exp(1j*k)*H1 - exp(-1j*k)*H1.H)
    
    return(Htot,Hcurrent)

def rajzLandau_E(theta,dim,klim,Npoints):
    '''
    drawing the Landau spectrum
    '''
    # theta = 0.025
    # dim = 100

    # klim = 1.7
    # Npoints = 200

    kran=linspace(-klim,klim,Npoints)
    dat=[]
    for k in kran:
        dat.append(eigvalsh(Htot(k, theta, dim)[0]))

    dat = array(dat)   ###   ez FONTOS 
    figure(figsize=(xfig_meret,yfig_meret))

#  first 11 levels 
    #plot(kran,dat[:,0:11])    
    plot(kran,dat)    

    emin=eigvalsh(Htot(0, theta, dim)[0])[0]
    hbarom= eigvalsh(Htot(0, theta, dim)[0])[1]-eigvalsh(Htot(0, theta, dim)[0])[0]
    
    en = [emin, emin+hbarom,emin+2*hbarom, emin+3*hbarom]
    labels = [r'$\frac{1}{2}\,\hbar \omega$', r'$\frac{3}{2}\,\hbar \omega$',
          r'$\frac{5}{2}\,\hbar \omega$', r'$\frac{7}{2}\,\hbar \omega$']
    yticks(en, labels , rotation='horizontal',fontsize=legend_meret)

    xlabel(r'wavenumber $k$',fontsize=xylabel_meret)
    ylabel(r'$E(k)$',fontsize=xylabel_meret); 

    xlim(-klim,klim)
    emax=eigvalsh(Htot(0, theta, dim)[0])[11]
    ylim(eigvalsh(Htot(0, 0, dim)[0])[0],-3.5)

    cim = 'Landau levels for strip' + ', $\Theta$ =' +str(theta) + ', W = '+str(dim)
    title(cim,fontsize=legend_meret)

 
    grid();
    
    #return();

def rajz_current(theta,dim,kin,n1,n2):
    
    # theta = 0.025
    # dim = 100

    # kin = 1.2

    nlevel = [n1,n2]

    res = Htot(kin, theta, dim)
    Hop = res[0]
    H_curr = res[1]
    ertek,vektor=eigh(Hop) 

    curr1 = []
    curr2 = []
    for i in range(dim):
        curr1.append(vektor[i,nlevel[0]]*transpose(H_curr*vektor[:,nlevel[0]])[0,i])
        curr2.append(vektor[i,nlevel[1]]*transpose(H_curr*vektor[:,nlevel[1]])[0,i])
    
    #leg1= 'current density for n = '+ str(nlevel[0]) + ', k =' + str(kin)
    #leg2= 'current density for n = '+ str(nlevel[1]) + ', k =' + str(kin)
    leg1= 'n = '+ str(nlevel[0]) + ', k =' + str(kin)
    leg2= 'n = '+ str(nlevel[1]) + ', k =' + str(kin)

    figure(figsize=(xfig_meret,yfig_meret))

    plot(curr1,'ro-',label=leg1)
    plot(curr2,'bo--',label=leg2)

    kin = -kin
    
    res = Htot(kin, theta, dim)
    Hop = res[0]
    H_curr = res[1]
    ertek,vektor=eigh(Hop) 

    curr1 = []
    curr2 = []
    for i in range(dim):
        curr1.append(vektor[i,nlevel[0]]*transpose(H_curr*vektor[:,nlevel[0]])[0,i])
        curr2.append(vektor[i,nlevel[1]]*transpose(H_curr*vektor[:,nlevel[1]])[0,i])
    
    #leg1= 'current density for n = '+ str(nlevel[0]) + ', k =' + str(kin)
    #leg2= 'current density for n = '+ str(nlevel[1]) + ', k =' + str(kin)
    leg1= 'n = '+ str(nlevel[0]) + ', k =' + str(kin)
    leg2= 'n = '+ str(nlevel[1]) + ', k =' + str(kin)
 
    plot(curr1,'rs-',label=leg1)
    plot(curr2,'bs--',label=leg2)

    legend(loc='upper center',fontsize=legend_meret)

    xlabel(r'$y/a$',fontsize=xylabel_meret); 
    ylabel(r'${\hbar k|\psi_n(y)|}^2$',fontsize=xylabel_meret); 
    xlim(0,dim)
    
    cim = 'Current density for Landau levels' + ', $\Theta$ =' +str(theta) + ', W = '+str(dim)
    title(cim,fontsize=15)
    
    grid();

    # return();

Landau levels for a strip of square lattice

theta = 0.025
dim = 100

klim = 2.
Npoints = 200

rajzLandau_E(theta,dim,klim,Npoints);

interact(rajzLandau_E,
               theta=FloatSlider(min=0.,max=0.1,step=0.005,value=0.025,description='theta='),
               dim=IntSlider(min=20,max=200,step=10,value=100,description='W='),
               klim =FloatSlider(min=-pi,max=pi,step=0.1,value=1.2,description='k_lim='),
               Npoints=IntSlider(min=50,max=200,step=10,description='Npoints='));

Wave function at given wave number $k$

for the first few Landau levels (given by parameter 'nlevel')

theta = 0.025
dim = 100

kin = 1.2

ertek,vektor=eigh(Htot(kin, theta, dim)[0])

nlevel = [0,1]

psi1 = []
psi2 = []
for i in range(dim):
    psi1.append(abs(vektor[i,nlevel[0]])*abs(vektor[i,nlevel[0]])) 
    psi2.append(abs(vektor[i,nlevel[1]])*abs(vektor[i,nlevel[1]]))
    
leg1= 'wave fn. for n = '+ str(nlevel[0])
leg2= 'wave fn. for n = '+ str(nlevel[1])
    
figure(figsize=(xfig_meret,yfig_meret))

plot(psi1,'ro-',label=leg1)
plot(psi2,'bo--',label=leg2)
legend()

xlabel(r'$y/a$',fontsize=xylabel_meret); 
ylabel(r'${|\psi_n(y)|}^2$',fontsize=xylabel_meret); 


cim = 'Wave function for Landau levels' + ', $\Theta$ =' +str(theta) + \
', W = '+str(dim) + ', k =' + str(kin)
title(cim,fontsize=15)



grid();

theta = 0.025
dim = 100

kin = 0

ertek,vektor=eigh(Htot(kin, theta, dim)[0])

nlevel = [0,1]

psi1 = []
psi2 = []
for i in range(dim):
    psi1.append(abs(vektor[i,nlevel[0]])*abs(vektor[i,nlevel[0]])) 
    psi2.append(abs(vektor[i,nlevel[1]])*abs(vektor[i,nlevel[1]]))
    
leg1= 'wave fn. for n = '+ str(nlevel[0])
leg2= 'wave fn. for n = '+ str(nlevel[1])
    
figure(figsize=(xfig_meret,yfig_meret))

plot(psi1,'ro-',label=leg1)
plot(psi2,'bo--',label=leg2)
legend()

xlabel(r'$y/a$',fontsize=xylabel_meret); 
ylabel(r'${ |\psi_n(y)|}^2$',fontsize=xylabel_meret); 

cim = 'Wave function for Landau levels' + ', $\Theta$ =' +str(theta) + \
', W = '+str(dim) + ', k =' + str(kin)

title(cim,fontsize=15)
grid();

Current density across the strip

at a given wavenumber $k$

for the first few levels (given by parameters 'n1' and 'n2')

theta = 0.025
dim = 100

kin = 1.2
n1, n2= [0,1]

figure(figsize=(xfig_meret,yfig_meret));
rajz_current(theta,dim,kin,n1,n2);

/home/cserti/.local/lib/python3.4/site-packages/numpy/core/numeric.py:482: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

<matplotlib.figure.Figure at 0x7f3df4264588>

interact(rajz_current,
        theta=FloatSlider(min=0.,max=0.1,step=0.005,value=0.025,description='theta='),
        dim=IntSlider(min=20,max=200,step=10,value=100,description='W='),
        kin=FloatSlider(min=0,max=pi,step=0.1,value=1.2,description='kin='),
        n1=IntSlider(min=0,max=7,step=1,value=0,description='n1='),
        n2=IntSlider(min=0,max=7,step=1,value=1,description='n2='));

/home/cserti/.local/lib/python3.4/site-packages/numpy/core/numeric.py:482: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Spectrum and current density togehter:

def rajz_spectrum_current(theta,dim,kin,n):
    
    # theta = 0.025 
    # dim = 100

    # kin = 1.2

    # n ---> number of Landau nivo

    
    ######################################################################
    ###  plot Landau spectrum
    
    #figure(figsize=(xfig_meret,yfig_meret))
    
    klim = 2
    Npoints = 200

    kran=linspace(-klim,klim,Npoints)
    dat=[]
    for k in kran:
        dat.append(eigvalsh(Htot(k, theta, dim)[0]))

    dat = array(dat)   ###   ez FONTOS 
   
    figsize(11,5)
    figsize(14,7)
    
    subplot(1,2,1)  
    plot(kran,dat)    

    ertek = eigvalsh(Htot(kin, theta, dim)[0])
    plot(kin,ertek[n],'ro',markersize=10)
    plot(-kin,ertek[n],'bs',markersize=10)
    
    emin=eigvalsh(Htot(0, theta, dim)[0])[0]
    hbarom= eigvalsh(Htot(0, theta, dim)[0])[1]-eigvalsh(Htot(0, theta, dim)[0])[0]
    
    en = [emin, emin+hbarom,emin+2*hbarom, emin+3*hbarom]
    labels = [r'$\frac{1}{2}\,\hbar \omega$', r'$\frac{3}{2}\,\hbar \omega$',
          r'$\frac{5}{2}\,\hbar \omega$', r'$\frac{7}{2}\,\hbar \omega$']
    #yticks(en, labels , rotation='horizontal',fontsize=legend_meret)

    xlabel(r'wavenumber $k$',fontsize=xylabel_meret)
    ylabel(r'$E(k)$',fontsize=xylabel_meret); 

    xlim(-klim,klim)
    emax=eigvalsh(Htot(0, theta, dim)[0])[11]
    ylim(eigvalsh(Htot(0, 0, dim)[0])[0],-3.)
 
    cim1 = 'Energy dispersion'
    title(cim1,fontsize=20)
    
    grid();
    
    ######################################################################
    ###  plot current density

    subplot(1,2,2)
    
    
    res = Htot(kin, theta, dim)
    Hop = res[0]
    H_curr = res[1]
    ertek,vektor=eigh(Hop) 

    curr1 = []
    curr2 = []
    for i in range(dim):
        curr1.append(vektor[i,nlevel[0]]*transpose(H_curr*vektor[:,nlevel[0]])[0,i])
        curr2.append(vektor[i,nlevel[1]]*transpose(H_curr*vektor[:,nlevel[1]])[0,i])
    
    #leg1= 'current density for n = '+ str(nlevel[0]) + ', k =' + str(kin)
    #leg2= 'current density for n = '+ str(nlevel[1]) + ', k =' + str(kin)
    leg1= 'n = '+ str(nlevel[0]) + ', k =' + str(kin)
    leg2= 'n = '+ str(nlevel[1]) + ', k =' + str(kin)

    
    plot(real(curr1),'ro-',label=leg1)
    plot(real(curr2),'bo--',label=leg2)

    res = Htot(-kin, theta, dim)
    Hop = res[0]
    H_curr = res[1]
    ertek,vektor=eigh(Hop) 

    curr1 = []
    curr2 = []
    for i in range(dim):
        curr1.append(vektor[i,nlevel[0]]*transpose(H_curr*vektor[:,nlevel[0]])[0,i])
        curr2.append(vektor[i,nlevel[1]]*transpose(H_curr*vektor[:,nlevel[1]])[0,i])
    
    #leg1= 'current density for n = '+ str(nlevel[0]) + ', k =' + str(kin)
    #leg2= 'current density for n = '+ str(nlevel[1]) + ', k =' + str(kin)
    leg1= 'n = '+ str(nlevel[0]) + ', k =' + str(kin)
    leg2= 'n = '+ str(nlevel[1]) + ', k =' + str(kin)
 
    plot(real(curr1),'rs-',label=leg1)
    plot(real(curr2),'bs--',label=leg2)
    
    
    legend(loc='upper center',fontsize=legend_meret)

    xlabel(r'$y/a$',fontsize=xylabel_meret); 
    ylabel(r'${\hbar k |\psi_n(y)|}^2$',fontsize=xylabel_meret); 
    
    ticklabel_format(style = 'sci', useOffset=False)
    
    xlim(0,dim)
    #ylim(-0.2,0.2)
    
    cim2 = 'Current density '
    title(cim2,fontsize=20)
  
  
    cimfo = 'Landau levels for strip' + '\n $B\sim \Theta $ =' +str(theta) + ', W = '+str(dim)
    
    grid();
    
    
    tight_layout();
  
    # shift subplots down:
    
    st = suptitle(cimfo,fontsize=20)  ###  a ket vizszintes abra kozos cime
    st.set_y(0.95)
    subplots_adjust(top=0.8)
    

theta = 0.025
dim = 100

kin = 1.7
n = 0

rajz_spectrum_current(theta,dim,kin,n);

#rajz_spectrum_current(theta,dim,kin,n)
    
interact(rajz_spectrum_current,
        theta=FloatSlider(min=0.,max=0.1,step=0.005,value=0.025,description='theta='),
        dim=IntSlider(min=20,max=200,step=10,value=100,description='W='),
        kin=FloatSlider(min=0,max=1.7,step=0.1,value=1.2,description='kin='),
        n=IntSlider(min=0,max=7,step=1,value=0,description='n='));